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Abstract 

Constitutive defense mechanisms are critical to the understanding of defense mechanisms in conifers because they 
constitute the first barrier to attacks by insect pests. In interior spruce, trees that are putatively resistant and susceptible to 
attacks by white pine weevil (Pissodes strobi) typically exhibit constitutive differences in traits such as resin duct size and 
number, bark thickness, and terpene content. To improve our knowledge of their genetic basis, we compared globally the 
constitutive expression levels of 17,825 genes between 20 putatively resistant and 20 putatively susceptible interior spruce 
trees from the British Columbia tree improvement program. We identified 54 upregulated and 137 downregulated genes in 
resistant phenotypes, relative to susceptible phenotypes, with a maximum fold change of 2.24 and 3.91, respectively. We 
found a puzzling increase of resistance by downregulated genes, as one would think that "procuring armaments" is the best 
defense. Also, although terpenes and phenolic compounds play an important role in conifer defense, we found few of these 
genes to be differentially expressed. We found 15 putative small heat-shock proteins (sHSP) and several other stress-related 
proteins to be downregulated in resistant trees. Downregulated putative sHSP belong to several sHSP classes and 
represented 58% of all tested putative sHSP. These proteins are well known to be involved in plant response to various kinds 
of abiotic stress; however, their role in constitutive resistance is not yet understood. The lack of correspondence between 
transcriptome profile clusters and phenotype classifications suggests that weevil resistance in spruce is a complex trait. 
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Introduction 

The white pine weevil (Pissodes strobi) is a major pest of 
North American forests (Drouin and Langor 1991; Alfaro 
1 994; Hamid et al. 1 995). The weevil primarily attacks Sitka 
spruce (Picea sitchensis), white spruce (Pi. glauca), and En- 
gelmann spruce (Pi. engelmannii), but it can also attack sev- 
eral other pine and spruce species and even Douglas fir 
(Pseudotsuga menziesii). Adults lay eggs in the bark below 
the terminal bud cluster, and larva feed on the terminal 
leader. Such attacks can lead to leader death and conse- 
quential stem deformation, which is an economic cost to 
the forest industry (Alfaro 1 994). Knowledge of the genetic 
mechanisms of weevil resistance in spruce would aid in de- 
veloping marker-assisted breeding strategies for spruce and 
add to our knowledge about the diversity of resistance 
mechanisms in the plant kingdom. 



In conifers as in other plants, resistance to insect pests 
involves both constitutive (pre-existing) and induced de- 
fenses. Constitutive defense mechanisms are both mechan- 
ical (resin ducts, parenchyma cells, and sclerenchyma) and 
chemical (oxalate crystals and accumulation of toxic or re- 
pellant molecules) (Hall et al. 201 1 ). Induced defenses form 
a second line of defense, operating during or after pest at- 
tack. They are generally more specific in their action and in- 
clude increases of resin flow and production of repellant or 
toxic chemicals or even de novo defenses (formation of trau- 
matic resin ducts, callus formation, synthesis of new chem- 
icals that are possibly specific to a given pest). Most workers 
regard induced defensive mechanisms to be the most impor- 
tant component of insect defense; however, constitutive re- 
sistance is less liable and easier to study and quantify in the 
context of quantitative genomics. 
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With regard to white spruce, several studies have identi- 
fied constitutive features of resistance. Resistant trees pos- 
sess a thinner bark, with a higher density of outer resin ducts 
and larger inner resin ducts (Tomlin and Borden 1994, 
1997b; Alfaro et al. 2004). In interior spruce (Pi. glauca- 
engelmannii complex), resistance is positively correlated 
with tree growth (both height and trunk diameter), 
although weevils prefer to oviposit in longer leader shoots 
(Kiss and Yanchuk 1991; King et al. 1997). Gerson and 
Kelsey (2002) analyzed piperidine alkaloids contents of 
resistant and susceptible families of Sitka spruce, but they 
did not find any correlation with resistance to weevil ovipo- 
siting. With regard to terpenoids, Nault et al (1 999) showed 
profiles to be good indicators of resistance in white spruce 
and Engelmann spruce. In Sitka spruce, resistant trees can 
show either a lower or a higher content of foliar terpenoids 
than susceptible trees, suggesting that they can use either 
repellency strategy (the tree try to repel the insects) or 
stealth strategy (the tree try to be less attractive to the 
insects; Tomlin etal. 1997). However, higher levels of a diter- 
pene (dehydroabietic acid) and two monoterpenes ((+)-3- 
carene and terpinolene) are associated with resistance in 
Sitka spruce (Robert et al. 2010). Following this study, Hall 
et al. (201 1) showed that the (+)-3-carene is produced by 
three different (+)-3-carene synthase genes. One was 
specific to resistant trees (PsTPS-3car2) f one was specific 
to susceptible trees (PsTPS-3car3), and one is expressed in 
both phenotypes (PsTPS-3car1). They concluded that 
(+)-3-carene is explained by the variation in gene copy 
number, in gene sequence, in protein expression levels, 
and in enzyme activity levels. 

The development of "omics" approaches and the devel- 
opment of several cDNA libraries within the Arborea I, II and 
Treenomix I, II spruce genome projects (http://www. arborea. 
ulaval.ca/; http://www.treenomix.ca; Pavy et al. 2005; 
Ralph, Yueh, et al. 2006; Ralph et al. 2008) opened insights 
into the nature of both constitutive and induced defense 
mechanisms in spruce. To date, most published studies have 
focused on induced defenses (Ralph, Yueh, et al. 2006; 
Lippert et al. 2007, 2009; Zulak et al. 2009; Robert et al. 
201 0; Hall et al. 201 1 ). These studies compare the biological 
response with various types of induction (methyl jasmonate 
and chitosan elicitation, white pine weevil and western 
spruce budworm herbivory, mechanical wounding) at the 
transcriptome, proteome, and/or metabolome levels. How- 
ever, induced and constitutive defenses are complementary 
and distinct defense mechanisms. Induced defenses take 
place when constitutive defenses have been defeated by 
an insect attack. Their primary function is to reinforce the 
constitutive defense mechanisms and add new barriers 
against the insect attack. Consequently, we might expect 
induced and constitutive defenses to have a different ge- 
netic basis. The purpose of this study was to investigate 
these differences. 



The comparison of resistant and susceptible trees at the 
global transcriptome level has not yet been conducted, and 
such a comparison can provide fundamental and perhaps 
unexpected findings about the basis of insect resistance 
in conifers. Here, we present a comparative study of gene 
expression in interior spruce (Pi. glauca-engelmannii com- 
plex) aimed to identify candidate genes involved in consti- 
tutive defense against white pine weevil. We used a set of 
180 trees previously ranked for resistance to this weevil by 
breeders in the British Columbia Ministry of Forests. Using 
a 17,825 member cDNA microarray, we compare gene ex- 
pression levels between the 20 most resistant trees and the 
20 most susceptible trees. Significantly upregulated and 
downregulated genes will identify a suite of genes involved 
in constitutive weevil resistance. Particular attention will be 
given to the putative small heat-shock proteins (sHSP) that 
evidently play an important role in constitutive defense. 



Materials and Methods 

Selection and Sampling 

As part of the British Columbia (BC) interior spruce tree 
breeding program (Experimental Project EP 670), 180 trees 
were selected in wild stands across the Prince George region 
of central BC (fig. 1). The parent tree selection criteria were 
largely height superiority, stem form, branch size, and crown 
shape. Their ages varied from 100 to 200 years. Open-pol- 
linated seeds were collected from each wild tree, and test 
seedlings for each parent tree were grown in nursery beds 
near Prince George. Progeny tests of all families were estab- 
lished in 1972 at Aleza Lake, near Prince George, and in 
1973 at three other sites: the Prince George Tree Improve- 
ment Station (PGTIS), Quesnel, and Barbie Lake. In the mid- 
1 980s, the PGTIS and Aleza Lake sites began to suffer severe 
and repeated attacks of white pine weevils. In 1988, 
presence or absence of weevil damage was recorded for 
all trees on both sites. Kiss and Yanchuk (1991) reported 
that family damage was consistent between the two sites 
(r = .71) and had a moderately strong genetic basis 

(^am,ly = 0 - 77 '^nd,v,duar°- 18 )- Kin 9 et aL < 1997 ) re P 0rted 

similar results in other BC interior spruce populations. Based 
on these results, it appears that parental resistant scores can 
be readily estimated from weevil damage on their proge- 
nies. In 2003, all families on both sites were ranked accord- 
ing to the number of damaged trees, and the observed 
damage was used to estimate resistance levels of the 180 
parent trees. In this study, the 20 least and 20 most dam- 
aged families were chosen as the resistant and susceptible 
families, respectively. 

In addition to collecting open-pollinated seed from the 
180 parent trees in the wild, scionwood (i.e., shoot tips) 
was collected from each tree, and all trees were cloned 
by grafting and established in clone banks at Vernon, Barnes 
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Fig. 1. — Parent trees' origin within the Prince George area. The color scale (S-R) indicates the level of resistance of the trees, from highly 
susceptible to highly resistant, blue to red, respectively Filled circles represent origin of the tree families used in the present microarray study. Open 
circles represent the origin of tree families not used in the microarray study but used for the resistance ranking (map layers from MapPlace Web site 
http://www.empr.gov.bc.ca/MINING/GEOSCIENCE/MAPPLACE/Pages/default.aspx). 



Creek (near Enderby, BC), and PGTIS. Samples used for ge- 
netic analysis in this study were collected from parent tree 
grafts at the Barnes Creek site. The use of cloned trees grow- 
ing in the same location instead of wild trees located across 
a vast geographic area removes bias due to different envi- 
ronmental growth conditions. 

RNA Extraction and Microarray Profiling 

Bark samples were collected from lateral shoots of the trees 
from the Barnes Creek clone bank. Total RNA was extracted 
following Kolosova et al. (2004). RNA quantity and quality 
were assessed by measuring spectral absorbance between 
200 and 350 nm and by visual assessment on a 1 % agarose 
gel. cDNA synthesis was completed for each sample indepen- 
dently using Superscript II Reverse Transcriptase (Invitrogen) 
with an oligo dT1 2-1 8 primer. cDNA samples were hybrid- 
ized using 3DNA Array 350 Expression Array Detection Kit 
(Genisphere) onto the Treenomix Spruce cDNA microarray 
(21. 8K version) comprising 18,725 unique elements. A bal- 
anced design with dye swaps was used to make direct com- 
parison of gene expression levels of resistant and susceptible 
trees. Each resistant tree was randomly contrasted with a sus- 
ceptible tree. 



Statistical Analysis 

Slides were scanned, and spot intensity was quantified using 
ImaGene 6.0.1 software (BioDiscovery, Inc., El Segundo, 
CA). To correct for background intensity, the lowest 10% 
of median foreground intensities per subgrid was subtracted 
from the median foreground intensities. Data were then 
normalized slide by slide, by variance stabilizing normaliza- 
tion to compensate for nonlinearity of intensity distributions 
(Huber et al. 2002). A linear mixed effects model was fit to 
the data taking account of both resistance/susceptibility and 
dye effects. Fold change (FC) and Pand 0 values were com- 
puted for all genes. Genes were considered to have a signif- 
icant differential expression level when their Pvalue is below 
0.05 and their fold change above 1 .5. 

Heat map and cluster analysis were performed on genes 
with P < 0.05 and FC > 1.5. Individuals and genes were 
clustered with Pearson correlation index and Spearman cor- 
relation index, respectively. Dendograms were drawn using 
the "hclust" function in R Script. 

To identify major themes appearing among the differ- 
entially expressed genes, we used the software Blast2Go 
(Conesa et al. 2005; Gotz et al. 2008) to test for statistical 
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overrepresentation of gene ontology terms (GO terms) 
among genes up- and downregulated. A more detailed func- 
tional categorization was performed using BlastXand tBIastX 
search versus viridiplantae database on National Center for 
Biotechnology Information (NCBI). We considered only re- 
sults with an E value lower than 10~ 10 . Given the number 
of differently expressed putative sHSP, a particular emphasis 
has been given to this protein family. tBIastN searches using 
protein sequences of known sHSP of Arabidopsis thaliana 
and Oryza sativa (Scharf et al. 2001; Siddique et al. 
2008; Sarkar et al. 2009) were performed over the whole 
microarray to identify putative members of the sHSP family. 
Sixty-one representative sequences of the 16 known sHSP 
classes from Ar. thaliana (Scharf et al. 2001; Siddique 
et al. 2008), Populus trichocarpa (Waters et al. 2008), 
and O. sativa (Sarkar et al. 2009) were added to this sequen- 
ces data set. Sequences were first aligned using the online 
version of PROMALS3D (Pei et al. 2008) and then optimized 
manually. The evolutionary distances were computed using 
the Poisson correction method (Zuckerkandl and Pauling 
1965). Phylogenetic relationships were inferred based on 
amino acid sequences using the Neighbor-Joining method 
to determine the exact class of each sequence. Only the con- 
served C-terminal sequences have been considered (see 
supplementary tables 1 and 2, Supplementary Material on- 
line). The reliability of the inferred tree was tested by boot- 
strap analysis with 1,000 replicates (Felsenstein 1985). 

Raw data and normalized data are uploaded to the 
Gene Expression Omnibus with accession number 
GSE27476 (http://www.ncbi.nlm.nih.gov/geo/query/acc. 
cgi?acc=GSE27476). Sequences for array clones are found 
in NCBI using the clone IDs given in Table 2 and Table 3 
and supplementary table S1 (Supplementary Material online). 

Results 

Resistance Levels 

The percentage of trees damaged by weevils was signifi- 
cantly higher among susceptible trees (68%) than among 
resistant trees (2 1 %; P< 2.2 x 10~ 16 ; fig. 2). No difference 
was found between susceptible and resistant trees neither in 
size nor in survival. Supplementary table S1 (Supplementary 
Material online) summarizes the observed damages. These 
results show that we have a valid comparison of phenotypic 
differences between two classes of trees that differ in 
resistance to the weevil. 

Gene Expression Profiles 

Among the 18,725 genes on our microarray chip, 2,499 
showed a P value less than 0.05 for significant differences 
of gene expression between the two classes of trees that 
differ in resistance (table 1). The highest 0 value observed 
among these genes was 0.282 but only one gene showed 



% of damaged trees 




Rank 



Fig. 2. — Percentage of damaged trees among progenies. Proge- 
nies are ordered from the least damaged to the most damaged. 
Resistant and susceptible families are located on the left and on the 
right, respectively. White bars and black bars show selected families for 
the present study. 

a 0 value less than 0.05. FC were low with the maximums 
FC of 2.24 and 3.91 in upregulated and downregulated 
genes, respectively (table 1 and fig. 3). Consequently, we 
considered gene expression to be significantly different if 
the P value was less than 0.05 and FC was greater than 
1.5. With such a rigorous criteria, we identified 54 genes 
as upregulated and 137 genes as downregulated, in resis- 
tant trees compared with susceptible trees, for a total of 
191 significant genes. 

As a further verification of differential gene expression, 
we performed cluster analysis and heat map based on 
the 191 significant genes (fig. 4). The cluster analysis indi- 
cates two groups, however, they do not match the resistant/ 
susceptible classification; cluster #1 contained 1 1 suscepti- 
ble trees, whereas cluster #2 contained 9 susceptible trees 
and 20 resistant trees. There is no evidence of a link between 
the resistance levels and the classification of susceptible 
trees in two distinct groups. The heat map (fig. 4) confirms 
the differences in gene expression profiles between the two 
clusters and suggests no difference between susceptible and 
resistant trees in cluster #2. Genes cluster in two main 
groups: 1) downregulated genes and 2) upregulated genes. 

To find differences that might exist between resistant 
trees and susceptible trees of clusters #1 and #2 (fig. 4), 
we performed a complementary analysis. We fitted the data 
as previously described to a mixed linear model but consid- 
ered three groups of trees: group S1 = cluster #1 (S-1 57- 
162, S-154-135, S-163-166, S-160-176, S-164-163, 
S-1 65-65, S-162-111, S-174-128, S-1 59-43, S-1 55-62, 
S-1 69-72), group S2 = susceptible individuals of cluster 
#2 (S-1 70-1 07, S-1 76-1 33, S-161-60, S-1 56-1 03, S-1 58- 
131, S-1 67-95, S-1 73-1 17, S-1 79-1 05, S-1 66-1 30, see 
fig. 4), and group R = resistant individuals (of cluster #2). 
This approach is not compatible with our experimental de- 
sign as this analysis consists of three groups, and the exper- 
imental design was made to compare two groups. Hence, 
individuals are not properly balanced over dyes and groups. 
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Table 1 

Summary of f-Test Comparisons between Resistant and Susceptible 
Trees (= reference) 



18,725 Analyzed 




Down- 


Genes 


Upregulated 


regulated 


Genes with P value < 0.05 


1,225 


1,274 




(FDR = 28.2%) 


(FDR = 28.2%) 


Genes with FC > 1 .5 


60 


151 


Maximum FC 


2.24 


3.91 


Significant genes 


54 


137 


(P< 0.05 and FC > 1.5) 







Note. — FDR, false discovery rate. 



Moreover, this statistical approach is not adequate as we 
predefine groups according to their gene expression profiles 
prior to the statistical comparisons based on the gene ex- 
pression profiles. So results should be taken with caution. 
Only 30 genes are significantly differently expressed (FC 
up to 3.52) between group R and group S2 according to 
the criteria P < 0.05 and FC > 1.5 but with a 0 value of 
1 (table 2). This tends to confirm the low levels of difference 
between these groups. By contrast, the observed differen- 
ces between group S1 and group R are high with 274 up- 
regulated (FC up to 10.05) and 430 downregulated genes 
(FC up to 3.40) in group S1. 

Functional Characterization 

Using Blast2Go, we tested the occurrence of overrepre- 
sented GO terms among the set of significant genes arising 
from the comparison of resistant and susceptible trees com- 
pared with the entire microarray. Among the biological pro- 
cesses, only a few categories were overrepresented (fig. 5): 
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Fig. 3. — Smoothed densities color representation of volcano plot, 
showing the differential expression levels of 18,725 genes between 
resistant and susceptible trees. Significant downregulated and upregu- 
lated genes are shown in blue and red, respectively. 



"response to hydrogen peroxide," "response to heat," and 
"response to high light intensity" and several higher catego- 
ries. All belong to the wider category "response to stimulus." 
Among cellular components, the only overrepresented cate- 
gory is "microtubule-associated complex." Among molecular 
functions, the two lowest overrepresented categories are 
"Rho guanyl-nucleotide exchange factor" and "microtubule 
motor." Although the trees were not stimulated, the over- 
represented GO terms suggest that differentially expressed 
genes are involved in stress or stimulus responses, but their 
molecular functions remain obscure. 

To complete analysis of the GO terms, BlastX and tBIastX 
searches were preformed against Viridiplantae on NCBI to 
deduce the functions of these putative genes, using E values 
less than 1 0- 1 0. One hundred and six clones gave no results 
or matched sequences with unknown functions. We did find 
85 matches with annotations using either BlastX or tBIastX. 
Genes with significant blast results are presented in Table 3. 
Differentially expressed genes belong to various gene families 
with few apparent links, except for putatively stress-related 
genes (including the putative sHSP). Three genes were 
annotated as putative transcription factors and three genes 
are annotated as part of putative transposable elements, but 
their possible function here is unknown. 

Of the 191 genes either up- or downregulated between 
resistant and susceptible trees, we found very few differen- 
tially expressed genes to be putatively involved in phenylpro- 
panoid and terpenoid metabolisms. Only four genes were 
putatively assigned to the terpenoid metabolism: one puta- 
tive cytochrome P450, two putative delta-selinene— like 
synthases that were downregulated, and one putative zeatin 
O-glucosyltransferase that was upregulated. Eight to nine genes 
were putatively directly related to phenylpropanoid metabolism: 
a putative UDP-glycosyltransferase, a putative laccase, two 
putative phenylcoumaran benzylic esther reductase, a puta- 
tive zeatin O-glucosyltransferase, a putative caffeic acid 
O-methyltransferase, a putative Flavonol 4'-sulfotransferase 
and a putative cytochrome P450, and eventually the putative 
transcription factor (MYB16) that might be linked to phenyl- 
propanoid or terpenoid metabolism (Bedon et al. 2007). 

Differential Expression of sHSP and Stress-Related 
Proteins 

Of the 26 putative sHSP printed on our microarray chips, 1 5 
were downregulated in resistant trees. We compared their 
sequences with An thaliana, Po. tnichocanpa, and Zea mays 
sHSP sequences allowing class determination of the majority 
of these genes (fig. 6). The phylogeny is congruent with pre- 
vious classifications of sHSP (Scharf et al. 2001; Siddique 
et al. 2008; Waters et al. 2008; Sarkar et al. 2009) with 
the exception of Os21.8ER, which was previously character- 
ized as a member of the endoplasmic reticulum group of 
sHSP but clustered here with Os 18.8 of the cytoplasmic 
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Fig. 4. — Heat map of the 191 significantly differently expressed genes between susceptible and resistant trees to the white pine weevil. Blue and 
red squares at the top of the heat map indicate susceptible and resistant trees, respectively. Tree labels are indicated at the bottom as follow: the tree 
phenotype (R = resistant, S = susceptible), the family rank in progeny tests for resistance (1 = the most resistant; 1 79 = the most susceptible), and then 
the family number. 



class X. Most of these spruce sHSP sequences cluster within 
the classes of sHSP previously identified in Ar. thaliana, Po. 
trichocarpa, and Z mays. Those that failed to cluster might 
belong to new sHSP classes. 

As in other species, the most diverse class of putative 
sHSP in spruce is the nucleocytoplasmic class I, represented 
by seven putative clones (WS0052_F03, WS00923_A06, 
WS0061JM21, WS0262_N22, IS0014J.07, WS0261_O21, 
and WS00823_I.11; fig. 6). Nucleocytoplasmic classes II 
and III are represented by two putative clones 
(WS0266_N22 and WS00825_O14) and one putative clone 
(WS00815_E02), respectively. WS0058_F08 putatively be- 
longs to the peroxisomal class and WS0063_C15 and 
WS00919J02 putatively belong to the endoplasmic retic- 
ulum class. WS0087J23, WS0058_B04, and 
WS00925_H13 do not cluster within any classes of either 
reference species. They may belong to a new class, specific 



to conifers. Six clones are found within a clade consisting of 
mitochondrial (group I) and chloroplastic sHSP. 
IS0014_C09 and WS0263_F23 unambiguously cluster 
within the mitochondrial group I of sHSP. Similarly, 
WS0063_G17 and WS00924_D21 unambiguously cluster 
within chloroplastic sHSP. Because WS0064_K01 and 
WS0061_H08 are branched between mitochondrial group 
I and chloroplastic sHSP within the large clade consisting of 
both mitochondrial and chloroplastic sHSP, they cannot be 
assigned with high confidence to either class. 
WS0092_E18, WS00826_O04, and WS0054_N08 do not 
match any known class of sHSP. Nevertheless, they are pu- 
tatively related to the cytosolic classes V, VI, and VII, respec- 
tively, and are tentatively assigned to these groups of sHSP. 
WS00930_B15 cannot be assigned to any sHSP class be- 
cause the clone sequence is too short, even though tBIastN 
and tBIastX searches place it as a putative sHSP. 
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Table 2 

Summary of f-Test Comparisons between Resistant and Susceptible Trees of Clusters #1 and #2 (= references) 





Resistant (20) vs. Group S1 (11) 


Resistant (20) vs. Group S2 (9) 


18,725 Analyzed Genes 


Upregulated 


Downregulated 


Upregulated 


Downregulated 


Genes with P value < 0.05 

Genes with FC > 1 .5 
Maximum FC 
Significant genes 
(P< 0.05 and FC > 1.5) 


1,778 

(FDR = 18.9%) 
326 

3.39 
274 


1,709 

(FDR = 18.9%) 
482 

10.04 
430 


305 

(FDR = 100%) 
79 
2.84 
15 


337 

(FDR = 100%) 
56 

3.22 
15 



Note. — FDR, false discovery rate. 



The downregulated putative sHSP belong to several classes 
working in different cellular compartments: nucleocytoplas- 
mic (nine putative sHSP of classes I — VI), endoplasmic reticu- 
lum (two putative sHSP), peroxisome (one putative sHSP), and 
chloroplast (one putative sHSP). Two of the downregulated 
putative sHSP could not be assigned to a particular class 
and operate in an unknown cellular compartment and seem 
to belong to the new sHSP class. In addition to these putative 
sHSP, 1 4 putative stress-related proteins of various gene fam- 
ilies are differentially expressed (1 2 downregulated and 2 up- 
regulated in resistant trees), including 3 putative heat-shock 
proteins and at least 2 putative universal stress proteins. 



Discussion 

Differences between Resistant and Susceptible 
Trees 

Our comparison gene expression for 1 8,725 genes between 
20 susceptible and 20 resistant trees found 54 upregulated 
genes and 137 downregulated genes in resistant trees as 
compared with the susceptible trees. As presented in the 
introduction, several studies have shown that differences 
exist between resistant and susceptible phenotypes at the 
morphological, chemical, and genetic levels. Moreover, pre- 
vious studies have shown several hundred genes are 



Biological process 




/ 








response to stimulus 
P value: 1.6 x 10 5 






response to abiotic 
stimulus 
P value: 2.0 x 10 7 



response to 
stress 
P value: 1.5 x 10 6 



response to oxidative 
stress 
P value: 4.8 x 10 5 



response to 
temperature 
stimulus 
P value: 7.7 x 10 1 



response to heat 
P value: 2.4 x 10 11 



response to reactive 
oxygen species 
P value: 3.3 x 10 7 



31 



3~ 



response to 
radiation 
P value: 7.5 x 10 4 



response to hydrogen 
peroxide 
P value: 3.6 x 10 6 



3~ 



response to light 

stimulus 
P value: 7.5 x 10 4 



response to light 

intensity 
P value: 3.0 x 10 6 



response to high light 
intensity 
P value: 2.1 x 10' 7 




GTPase regulator 

activity 
P value: 5.3 x 10" 4 



motor activity 
P value: 2.4 x 10 4 




small GTPase 
regulator activity 
P value: 5.3 x 10 4 




guanyl-nucleotide 
exchange factor 

activity 
P value: 2.9 x 10 6 



microtubule motor 

activity 
P value: 3.9 x 10 5 



Ras guanyl-nucleotide 
exchange factor 

Activity 
P value: 7.2 x 10" 7 



Rho guanyl-nucleotide 
exchange factor 

Activity 
P value: 7.2 x 10 7 



Cellular component 




/ 








microtubule associated 






complex 






P value: 


2.5 x 10' 5 





Fig. 5. — Significantly overrepresented GO terms of genes among significant upregulated or downregulated genes between susceptible and 
resistant trees, respectively. Fisher's exact tests with multiple testing corrections were performed using Blast2GO software. Only GO categories with false 
detection rate lower than 0.05 are shown. 
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Table 3 

Functional Categorization of Differentially Expressed Genes (P value < 0.05 and FC > 1.5) 



gi # Clone ID BlastX E Value tBIastX £ Value FC P Value Q Value 



Cytochrome P450 family and Terpenoid metabolism 



gi|49040156 


WS0101_H07 


Abscisic acid 8' -hydroxylase 


6 


X 


10" 




Cytochrome P450, putative [Ricinus communis] 


8 


X 


10 




1.6 


0.006 


0.144 






cytochrome P450 [Lactuca sativa] 




















1.55 






gi|49062217 


WS00712_A10 


Delta-selinene-like synthase 
[Picea sitchensis] 


6 


X 


10 


-70 


Delta-selinene-like synthase [Pi. sitchensis] 


4 


X 


10 


-101 


0.024 


0.221 


gi|69352521 


WS00929_B22 


Delta-selinene-like synthase 


2 


X 


10 


-119 


Delta-selinene-like synthase [Pi. sitchensis] 


7 


X 


10 


-162 


1.96 


0.039 


0.258 






[Pi. sitchensis] 


























Phenylpropanoid metabolism 


























gi|49057427 


WS0263_L06 


Caffeic acid O-methyltransferase 


2 


X 


10" 


-56 


O-methyltransferase [R. communis] 


8 


X 


10" 


-61 


1.76 


0.007 


0.154 






[Pinus pinaster] 




















1.75 






gi|49014799 


IS0012_L15 


Cytochrome P450— like protein 
[Arabidopsis 

thr)Hr)rir)] 
LI IOIIOI /C7J 


5 


X 


10" 


-23 


Cytochrome P450 [Populus trichocarpa] 


8 


X 


10" 


-24 


0.001 


0.102 


gi|49059256 


WS0071_C13 


UDP-glucosyltransferase 3 
[Pueraria montana var. lobata] 


1 


X 


10" 


-25 


Glucosyltransferase-8 [Vigna angularis] 


8 


X 


10" 




1.72 


0.034 


0.248 


gi|70636503 


WS00730_B15 


Laccase [Pinus taeda] 


4 


X 


10" 




Laccase (Lac7) [Pin. taeda] 


1 


X 


10" 






0.030 


0.240 


gi|49043266 


WS01011J14 


Phenylcoumaran benzylic ether reductase 
homolog TH6 [Tsuga heterophylla] 


2 


X 


10" 




Phenylcoumaran benzylic ether reductase 
homolog TH6 mRNA [T. heterophylla] 


2 


X 


10" 


-15 


:s 


0.004 


0.126 


gi|49056257 


WS0058_F16 


Phenylcoumaran benzylic ether 


2 


X 


10" 




Phenylcoumaran benzylic ether reductase- 


4 


X 


10" 






0.016 


0.199 






reductase-like protein [Po. trichocarpa] 










like protein [Po. trichocarpa] 










1 .56 






gi|49025769 


WS00928_F16 


Steroid sulfotransferase 1 [Brassica napus] 


3 


X 


10" 




Flavonol 4' -sulfotransferase, putative [R. 


2 


X 


10" 




0.002 


0.113 
















communis] 
















sHSP 




























gi|49056795 


WS0261_O21 


Chaperone [Agave tequilana] 


1 


X 


10" 


-49 


Heat-shock protein 18.2 [Ar. thaliana] 


9 


X 


10 


-56 


2.64 
3.91 


0.005 


0.140 


gi|49138681 


WS00823_L1 1 


Small heat-shock protein [Pseudotsuga 


1 


X 


10" 


-69 


Low-molecular weight heat-shock protein 


9 


X 


10" 


-92 


0.000 


0.096 






menziesii] 










[Ps. menziesii] 










2.01 






gi|4901 5440 


IS0014_L07 


Heat-shock protein 17.5 [Malus x 


1 


X 


10 


-54 


Small heat-shock protein [Malus x domestica] 


2 


X 


10 


-62 


0.008 


0.158 






domestica] 




















2.35 






gi|49054183 


WS0052_F03 


Small heat-shock protein [Ps. menziesii] 


1 


X 


10" 


-18 


HSP18.2 gene for 18.2 kDa heat-shock protein 


4 


X 


10 


-54 


0.009 


0.164 
















[Ar. thaliana] 










2.25 
1.7 

1 7 
I . / 






gi|49024180 


WS00923_A06 


Small heat-shock protein [Ps. menziesii] 


3 


X 


10" 


-54 


Heat-shock protein 18.2 [Ar. thaliana] 


8 


X 


10 


-60 


0.012 


0.182 


gi|49140326 


WS00825_O14 


Heat-shock protein 17.0 [Picea glauca] 


2 


X 


10" 


-12 


Heat-shock protein 17.0 [Pi. glauca] 


8 


X 


10" 


-78 


0.008 


0.161 


gi|49131870 


WS00815_E02 


Small heat-stress protein class CHI 


5 


X 


10" 


-32 


17.5 kDa class II heat-shock protein mRNA 


2 


X 


10" 


-24 


0.007 


0.157 






[Lycopersicon peruvianum] 










[Zea mays] 










1.5 
2.51 






gi|49023311 


WS00919J02 


Heat-shock protein, putative [R. communis] 


7 


X 


10" 


-37 


Heat-shock protein, putative [R. communis] 


4 


X 


10" 


-39 


0.001 


0.096 


gi|49016363 


WS0063_C15 


Heat-shock protein, putative [R. communis] 


1 


X 


10 


-29 


Small heat-shock protein (hsp21.4) mRNA 


4 


X 


10 


-31 


0.004 


0.123 
















[Cyclamen persicum] 










,,e 






gi|49056249 


WS0058_F08 


Peroxisomal small heat-shock protein 


3 


X 


10" 


-30 


Cytosolic class I small heat-shock protein 


4 


X 


10 


-30 


0.002 


0.111 






[Glycine max] 










HSP17.5 (hspl 7.5 gene) [Castanea sativa] 










1.93 






gi|49016448 


WS0063_G17 


Heat-shock protein [Ammopiptanthus 
mongolicus] 


6 


X 


10" 


-54 


Heat-shock protein (hsp) [Am. mongolicus] 


4 


X 


10" 


-59 


0.012 


0.183 



Table 3 Continued 



gi# 



Clone ID 



BlastX 



E Value 



tBlastX 



£ Value FC P Value Q Value 



gi|49026225 


WS00930. 


_B15 


Class II cytoplasmic small-molecular 
weight heat-shock protein 17.1 
[Pi. glauca] 


4x10 


-09 


Class II cytoplasmic small-molecular weight 
heat-shock protein 17.1 (EMB29, SMW 
HSP17.1) mRNA [Pi. glauca] 


4 x 10~ 11 


3.63 


0.001 


0.102 


gi|70654578 


WS00826. 


_O04 


17.3 kDa class I heat-shock protein 


1 x 10- 


-1 1 


No match 




„ 


0.001 


0.102 








[Gl. max] 










2.02 






gi|49056161 


WS0058_B04 


Hsp20.1 protein [Sola n u m peruvian urn] 


6 x 10~ 


-40 


Small heat-shock protein of cytosolic class I 


1 x 10- 44 


0.006 


0.145 














[Funaria hygrometrica] 




1.78 






gi|49025288 


WS00925. 


_H13 


17.7 kDa heat-shock protein [Helianthus annuus] 


1x10" 


-32 


17.7 kDa heat-shock protein gene [H. annuus] 


3 x 10- 32 


0.034 


0.248 



Heat-shock proteins family and other stress-related proteins 

gi|49047895 WS01 027_A08 Usp: Universal stress protein-like protein 4x10 

[Astragalus sinicus] 

gi|49021075 WS00914_D23 USP-like protein [As. sinicus] 9x 
gi|49025415 WS00926_B20 Stress-induced protein sti 1 -like protein 6x 

[Ar. thaliana] 

gi|49057553 WS0264_A17 Heat-shock 70 kD protein [Gl. max] 4x10 
gi|49015272 IS0014_D03 No match 



10 - z/ 

1Q -86 



gi|49023568 WS00920_E06 Chaperonin CPN60-2, mitochondrial (HSP60-2) 2 x 

[Cucurbita cv. Kurokawa Amakuri] 

gi|49051850 WS0018_D18 Hypothetic chloroplast chaperonin 21 [Vitis vinifera] 4x 

gi|49056969 WS0262_G19 Metallothionein-like protein [Sesamum indicum] 7x 

gi|4901 7024 WS0091J15 Ethylene-responsive protein, putative [Ar. thaliana] 1 x 

gi|4901 8081 WS0091 2J05 Jasmonate ZIM-domain protein 1 [So. lycopersicum] 2 x 

gi|70621 372 WS02610_H19 Arsenite-inducible RNA-associated protein aip-1, 4x 

putative [R. communis] 

gi|68771533 WS0064_H09 Alcohol dehydrogenase [Pinus banksiana] 1x 
gi|49023269 WS00919_G06 ATERDJ3A; oxidoreductase; putative DnaJ 1x 

protein [Ar. thaliana] 

gi|49132946 WS00816J13 Hypothetical water stress-induced protein 5x 

[Ps. menziesii] 

Transposable elements 



-82 
-15 



10 -io 

10 
10 

10~ 15 

10~ 12 

10- 47 

10 
10 

10 39 



-66 
-31 



ER6 protein (ethylene-inducible) [Solanum 

lycopersicum] 
USP-like protein mRNA [As. sinicus] 
Stress-inducible protein, putative [Ar. thaliana] 

Heat-shock protein [R. communis] 

Hsp90 mRNA for heat-shock protein 90 [Oryza 

sativa Japonica Group] 
Gland development-related protein 19-like mRNA 

[Gossypium hirsutum] 
cp1 0-like protein (CLP) mRNA [Gh. hirsutum] 
Seed-specific metallothionein-like protein (MT) 

gene [Se. indicum] 
Universal-stress protein (USP) family protein 

[Ar. thaliana] 
Jasmonate ZIM-domain protein 1 mRNA [So. 

lycopersicum] 
Arsenite-inducible RNA-associated protein 

aip-1, putative [R. communis] 
Alcohol dehydrogenase [Pin. banksiana] 
Heat-shock protein-binding protein, putative 

[R. communis] 
Galactinol synthase 1 [Po. trichocarpa x Populus 

deltoides] 



4 x 10 

3 x 10" 36 
1 x 10~ 102 

5 x 10" 54 

3 x 10~ 36 

1 x 10" 15 

4 x 10" 92 

2 x 10~ 15 

4 x 10~ 19 

1 x 10~ 10 

4 x 10 48 

2 x 10 

5 x 10 

2 x 10 39 



0.003 0.117 



1.8 
1.58 

2.39 
1.76 

0.58 

0.57 

:; 

0.63 
0.49 



0.006 
0.001 

0.001 
0.011 



-86 
-35 



1.59 
0.58 

3.31 



0.001 
0.027 

0.006 

0.023 

0.001 

0.003 
0.005 



0.148 
0.102 

0.102 
0.175 



0.003 0.120 



0.100 
0.230 

0.148 

0.220 

0.100 

0.118 
0.130 



0.003 0.117 



gi|491 36084 


WS00820_G23 


Copia-like retrotransposable element [Ar. thaliana] 


2 x 10~ 40 


Genes for S-locus F-Box protein c, Sc-Rnase 


1 x 10~ 43 




0.001 


0.106 










[Prunus dulcis] 




1.9 






gi|4901 8914 


WS0093_C16 


Copia-type polyprotein [Ar. thaliana] 


1 x 10~ 21 


Retrotransposon gtd1-12e3-re-5 [Glycine 


2 x 10~ 29 


0.015 


0.194 










tomentella] 




1.6 






gi|4901 761 9 


WS00911_D13 


Integrase [Boechera divaricarpa] 


2 x 10- 18 


Retrotransposon PpRT6 RNaseH-like gene 
[Pinus pinaster] 


1 x 10~ 42 


0.002 


0.108 



Table 3 Continued 



gi# 



Clone ID 



BlastX 



E Value 



tBlastX 



£ Value FC P Value Q Value 



Transcription factors 
















gi|49123245 WS0032_G19 


No match 


R2R3-MYB transcription factor MYB16 


9 x 


10 55 


2.69 


0.001 


0.096 






[Pi. glauca] 






1.58 






gi|49016284 WS0062_O09 


MBF1C (MULTIPROTEIN BRIDGING FACTOR 1C); DNA 
binding/transcription coactivator/transcription factor 
[Ar. thaliana] 


3 x 1 (T 50 msh6-2 gene, exon 1 to 1 7 [Ar. thaliana] 


2 x 


1Q -56 


0.017 


0.203 


gi|49015214 IS0014_A07 


Transcription initiation factor iib, putative [R. communis] 


5 x 10~ 13 Transcription initiation factor iib, putative 
[R. communis] 


3 x 


10" 15 




0.002 


0.112 



Other 

gi|49059326 WS0071_G04 

gi|49021 738 WS00915_F01 

gi|49025500 WS00926_F17 

gi|49042416 WS0108_M06 



gi|49055173 

gi|49022557 

gi|49052167 

gi|49025135 

gi|49017248 

gi|49025673 

gi|49040869 

gi|69354546 

gi|70634833 
gi|49045682 

gi|49018587 
gi|49019794 

gi|69359064 

gi|49023662 



WS0055_D11 

WS00917_H17 

WS00110_D02 

WS00922_M07 

WS00910_B07 

WS00928_B05 

WS0104J05 

WS00933_K09 

WS00724_G03 
WS01018_L23 

WS0092_C13 
WS0097_C17 

WS00937_N04 

WS00920J02 



AAA+-type ATPase (ISS) [Ostreococcus tauri] 
Alpha-glucan phosphorylase [Ar. thaliana] 

AT3G07090 [Ar. thaliana] 

ATCNGC4 (CYCLIC NUCLEOTIDE-GATED CATION 

CHANNEL 4); calmodulin 
binding/cation channel/cation transmembrane transporter/ 

cyclic nucleotide binding [Ar. thaliana] 
ATPP2-A4 (Phloem protein 2-A4); carbohydrate binding 

[Ar. thaliana] 

ATRBL14 (ARABIDOPSIS RHOMBOID-LIKE PROTEIN 14); 

zinc ion binding [Ar. thaliana] 
Cytoplasmic dynein light chain, putative [R. communis] 

Cytoplasmic dynein light chain, putative [R. communis] 

Cytoplasmic dynein light chain, putative [R. communis] 

Glycerophosphodiester phosphodiesterase [Z. mays] 

Hypothetical protein OsJ_1 431 5 [O. sativa Japonica 
Group] 

IQ calmodulin-binding region; Apoptosis regulator 

Bcl-2 protein, BAG [Medicago truncatula] 
Kinase, putative [R. communis] 
Late embryogenesis abundant protein [Pi. glauca] 

Metal ion-binding protein, putative [R. communis] 
Mitochondrial import inner membrane translocase 

subunit TIM 14 [Z. mays] 
Monovalent cation:proton antiporter, putative [R. 

communis] 

NADH:ubiquinone reductase subunit 2 [Beta 
vulgaris subsp. vulgaris] 



2 x 10 Zl 

5 x 10- 18 

1 x 10~ 31 

2 x 10~ 22 



8 x 10~ 11 

1 x 10~ 35 

2 x 10- 34 

2 x 10 29 

3 x 10~ 32 

4 x 10~ 61 

1 x 10~ 12 

2 x 10- 14 
3x10 



-40 
-49 



2x10 

2 x 10 41 
1 x 10~ 31 

1 x 10~ 15 

5 x 10- 18 



No match 

Alpha-1,4-glucan phosphorylase L isozyme 

[Cucurbita maxima] 
No match 

Putative ion channel, cngc4 [Ar. thaliana] 



No match 

ARABIDOPSIS RHOMBOID-LIKE PROTEIN 

14; ATRBL14 [Ar. thaliana] 
Cytoplasmic dynein light chain, putative 

[R. communis] 
Cytoplasmic dynein light chain, putative 

[R. communis] 
Cytoplasmic dynein light chain, putative 

[R. communis] 
Glycerophosphodiester phosphodiesterase 

[Z. mays] 
No match 

Bcl-2-associated athanogene-like protein 
[V. vinifera] 

Receptor-like kinase [Marchantia polymorpha] 
Late embryogenesis abundant protein (EMB6) 
[Pi. glauca] 

Metal ion-binding protein, putative [R. communis] 
Mitochondrial import inner membrane translocase 

subunit TIM 14 mRNA [Z mays] 
Monovalent cation:proton antiporter, putative 

[R. communis] 
NADH dehydrogenase subunit 2 [Cycas 

taitungensis] 



o.ooo 

3 x 10~ 41 157 0.012 



0.001 

1 x 10~ 20 0.012 



0.050 
0.180 

0.100 
0.180 




Table 3 Continued 



gi# Clone ID BlastX E Value tBIastX E Value FC P Value Q Value 



gi|49024651 


WS00924_F18 


Pirin, putative [Ar. thaliana] 


1 


X 


10" 


-12 


Pirin, putative [Ar. thaliana] 


1 X 


10 


-13 


0.61 
0.66 


0.002 


0.120 


gi|49016662 


WS0064_B23 


Protein phosphatase 2A regulatory A subunit 


2 


X 


10" 


-56 


Protein phosphatase 2A regulatory A subunit 


8 x 


io- 


-67 


0.028 


0.230 






[Lolium perenne] 
Putative hexose transporter [V. vinifera] 










mRNA [L. perenne] 
Sugar transporter, putative [R. communis] 














gi|49018463 


WS00913_L01 


1 


X 


10 


-17 


1 X 


10 






0.008 


0.160 


nii/iono c/i q 

gi|4yuDz_)4o 


WbUU I I l_rUo 


Putative neutral invertase [V. vinifera] 


z 


X 


1 rr 
I u 


-19 


nil gene for putative neutral invertase, exons 
1-4, clone 48C19 [V. vinifera] 


Z X 


I U 


-18 




r» mn 
U.UzU 


r» ~) 1 r» 
U.z I U 


gi|49016709 


\A/Qnn<^/i phi 

VV jUUD4_rU I 


Receptor-like protein kinase [Gl. max] 


7 
/ 


X 


1 rr 
I u 


-26 


Receptor-like protein kinase 1 [Gl. max] 


I X 


1 rr 
I U 


-34 


0.36 
0.56 
0.50 


U.UU I 


n 1 nn 
U. I uu 


gi|49017350 


WS00910_G03 


Retinol dehydrogenase 12 [Z. mays] 


9 


X 


10" 


-44 


Retinol dehydrogenase 12 [Z mays] 


1 X 


10 


-51 


0.020 


0.210 


gi|49024211 


WS00923_B17 


Retrotransposon protein [0. sativa (indica 


2 


X 


10 - 


-34 


Large subunit ribosomal RNA gene, partial 


2 x 


10 


-178 


0.009 


0.160 






cultivar-group)] 
RING-H2 finger protein ATL5A, putative 
[R. communis] 










sequence; chloroplast [Abies homolepis] 
RING-H2 finger protein ATL5A, putative 
[R. communis] 














g i|49045461 


WS01018_B20 


1 


X 


io- 


-21 


2 x 


10 


-24 




0.008 


0.160 


gi|4y i ozyby 
gi|49027191 
gi|4yu4ZDoy 


WbUUo I b_J I o 


Sor-like protein [Ginkgo biloba] 
Stem-specific protein TSJT1 [Z. mays] 


/ 


X 


1 rr 
I U 


-08 


Galactinol synthase [Coptis japonica] 
Stem-specific protein TSJT1 [Z mays] 


I X 


1 rr 
I U 


-40 


0.56 


U.UUU 


U. I UU 


WS00930_L23 


1 


X 


io- 


-26 


9 x 


io- 


-32 




0.005 


0.140 


\/\/cni no cc\q 

vvbu i uy_LUo 


Sterol regulatory element-binding protein site 2 


z 


X 


1 rr 
I U 


-72 


Sterol regulatory element-binding protein 


D X 


1 rr 
I U 


-82 


1.52 


r» nnc 
U.UUb 


U. I 4o 






protease, putative [R. communis] 










site 2 protease, putative [R. communis] 








2.21 






gi|4y 1 4U/ j i 


WbUUo I _N I U 


Transmembrane BAX inhibitor motif-containing 


7 


X 


1 0" 


-12 


Transmembrane BAX inhibitor 


2 x 


1 0~ 


-14 


0.002 


0. 1 09 






protein 4 [Z mays] 










motif-containing protein 4 [Z mays] 








1.57 






gi|4y i jo i jj 


VV jUUoZLMUy 


Transmembrane protein TPARL, putative [R. communis] 


A 

4 


X 


1 rr 
I U 


-45 


Transmembrane protein TPARL, putative 
[R. communis] 


Z X 


1 rr 
I U 


-76 


n r»r»3 

U.ULO 


U. I Z I 


gi|49025187 


WS00922_O17 


UBX domain-containing protein, putative [R. communis] 


9 


X 


io- 


-07 


UBX domain-containing protein, putative 
[R. communis] 


4 x 


10 


-20 




0.031 


0.240 


gi|49024645 


WS00924_F12 


Vacuole membrane protein, putative [R. communis] 


6 


X 


10 


-65 


Vacuole membrane protein, putative [R. 
communis] 


2 x 


10 


-71 


1.56 


0.004 


0.127 


gi|49044389 


WS01014_N23 


Zeatin O-glucosyltransferase [Phaseolus lunatus] 


1 


X 


io- 


-57 


Zeatin O-glucosyltransferase (ZOG1) [Ph. 
lunatus] 


1 X 


10 


-57 




0.006 


0.146 


gi|49143022 


WS0087_P05 


Zinc finger protein [Populus euphratica] 


2 


X 


io- 


-21 


NADH dehydrogenase [C. taitungensis] 


2 x 


10" 


-77 


0.65 
0.45 


0.008 


0.160 


gi|49024616 


WS00924_E06 


Zn-dependent hydrolases, including glyoxylases 


3 


X 


io- 


-29 


Metallo-beta-lactamase family protein [Ar. 


3 x 


io- 


-28 


0.000 


0.100 






[Z mays] 










thaliana] 








1.67 






gi|49055808 


WS0056_A24 


Putative callose synthase catalytic subunit [Go. 


1 


X 


io- 


-46 


Putative callose synthase catalytic subunit 


3 x 


10 


-52 


0.003 


0.121 






hirsutum] 










[Go. hirsutum] 








0.59 






gi|70630512 


WS0266_K23 


Dihydrofolate reductase-thymidylate synthase [Po. 
trichocarpa] 


7 


X 


10 


-31 


Difunctional dihydrofolate reductase-thymidylate 
synthase, putative [R. communis] 


7 x 


10 


-35 


0.002 


0.120 


gi|49020606 


WS0099J08 


Adipocyte plasma membrane-associated protein, 
putative [R. communis] 


7 


X 


io- 


-34 


Adipocyte plasma membrane-associated protein, 
putative [R. communis] 


1 X 


10 






0.001 


0.100 


gi|49063070 


WS0078_C13 


Germin-like protein [Ananas comosus] 


3 


X 


10- 


-41 


GLP5 (GERMIN-LIKE PROTEIN 5); manganese 
ion binding/nutrient reservoir [Ar. thaliana] 


8 x 


10- 


-68 




0.002 


0.110 



Note. — sHSP were first identified using BlastX and tBIastX search, and their class were further determined by phylogeny (see text). Colors green and red indicate downregulation and upregulation in resistant trees, respectively. gi#, Gl 
number of spruce clone in NCBI database. 
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Fig. 6. — Phylogenetic analysis of spruce sHSP. The tree was derived 
by Neighbor-Joining method with bootstrap analysis (1,000 replicates) 
from alignment of amino acid sequences of sHSP of rice, Arabidopsis, 
and poplar. Bootstrap values higher than 50% are shown next to the 



involved in induced defenses in both Sitka spruce and Nor- 
way spruce (Ralph, Yueh, et al. 2006; Lippert et al. 2009). 
Therefore, the number of differentially expressed genes (i.e., 
with FC higher than 1 .5) was expected to be greater than 
21 1 that found in this study (191 statistically significant). 

Such a low number of differentially expressed genes sug- 
gests that differences between resistant and susceptible 
phenotypes are linked more to variation in gene sequences, 
translation, and/or variation of catalytic efficiencies than to 
regulatory differences. Hall et al. (201 1) showed that differ- 
ences in (+)-3-carene levels can be explained by variation in 
1 ) the number of gene copies, 2) protein expression levels, 3) 
gene sequences, and 4) catalytic efficiencies. Such differen- 
ces can also be expected in other gene families, and the ob- 
served differences of gene expression levels may not explain 
all of the observed phenotypic differences. 

Another possible explanation for the low number of dif- 
ferentially expressed genes is that in conifers, several gene 
families are composed of a large number of closely related 
genes: terpenoid synthases (Martin etal. 2004; Keeling etal. 
2008), cytochrome P450 monooxygenases (Hamberger and 
Bohlmann 2006), dirigent proteins (Ralph, Park, et al. 2006; 
Ralph et al. 2007), and MYB transcription factors (Bedon 
et al. 2007, 2010). Therefore, we can expect that some 
spots of the microarray hybridize with transcripts of two, 
or even several, similar genes. In these cases, the observed 
gene expression levels are the average of the respective 
gene expression levels (i.e., upregulated genes cancel the 
effect of the downregulated genes). The low number of dif- 
ferentially expressed genes can also be linked to the exis- 
tence of disparate strategies of resistance (e.g., stealth or 
repellent) (see Phenotype Prediction and Efficiency of the 
Approach). 

Previous comparisons between resistant and susceptible 
trees have shown that resistant phenotypes in spruce are 
better "armed" to defend against weevils; however, these 
results are inconsistent. Tomlin and Borden (1994, 1997b) 
and Alfaro et al. (2004) found that resistant trees possessed 
more and larger resin ducts, whereas Tomlin et al. (1997) 
and Nault et al. (1999) reported no clear link between ter- 
pene profiles and resistance. Only one study suggested the 
existence of a stealth strategy (Tomlin et al. 1997). In the 
case that procuring "armaments" is the most common de- 
fense strategy, we might expect a majority of upregulated 
genes in resistant phenotypes. However, most of the differ- 
entially expressed genes in this study were downregulated 
(72%). This suggests that resistance could be linked more to 
a stealth strategy than to a repellent strategy. The silencing 
of certain genes may reduce the probability of detection and 
attack by weevils. Moreover, because resistance is useful 

branches. Phylogenetic analyses were conducted in MEGA4. EST clones 
ID of Picea are indicated in bold and underlined. Downregulated sHSP 
are indicated by a closed black circle. 
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only when weevils are present, the cost of a constant expres- 
sion of genes involved in resistance might be higher than the 
associated benefit. 

The comparison of resistant trees and the 1 1 susceptible 
trees of cluster #1 lead to a higher number of differentially 
expressed genes than the comparison of the 20 resistant 
and 20 susceptible trees. It suggests that more genes might 
show differences in constitutive expression levels. However, 
we cannot link the classification of the trees in three groups 
to a classification of phenotypes. Because this statistical ap- 
proach is not adequate, we will not talk more about these 
results and we just mention them as further analyses. 

Terpenoid and Phenylpropanoid Pathways: Few 
Genes Were Constitutively Differently Expressed 
in Resistant Spruces 

Only three differentially expressed genes have been found 
across the terpenoid metabolic pathways. Only two putative 
delta-selinene— like synthases are downregulated in resistant 
trees. In grand fir, delta-selinene synthase use farnesyl pyro- 
phosphate as substrate to produce 34 different sesquiter- 
pene olefins (Steele et al. 1998). The downregulated gene 
annotated as putative abscisic acid 8 '-hydroxylase belongs 
to the wide super family of cytochrome P450. This enzyme 
degrades abscisic acid into 8'-hydroxyabscic acid (Nambara 
and Marion-Poll 2005). Abscisic acid is an important terpe- 
noid phytohormone involved in many plant developmental 
processes and plant responses to environmental stress and 
pathogens (Seo and Koshiba 2002). In particular, abscisic 
acid regulates the opening of stomates and thus the loss 
of water in cells. Pei et al. (2000) showed abscisic acid also 
triggers an increase in cytosolic calcium in guard cells. In Pis- 
tia stratiotes, the Ca 2+ channels play an important role in 
calcium oxalate crystals formation (Volk et al. 2004). We 
might hypothesize that the reduced catabolism of abscisic 
acid is linked to an increase in the production of the toxic 
calcium oxalate crystals. However, more research is needed 
to confirm this hypothesis. 

There are seven differently expressed genes that can be 
putatively assigned to phenylpropanoid metabolism. First, 
a putative caffeic O-methyltransferase (COMT) is downre- 
gulated in resistant trees. This enzyme is known to be in- 
volved in methylation of precursors of both syringyl- and 
guaiacyl-lignin subunits in angiosperms (Baucher et al. 
2003; Do et al. 2007; Vanholme et al. 2008; Tu et al. 
2010). Several studies showed that downregulation of 
COMT leads to syringyl/guaiacyl-lignin ratio change or event 
suppression of syringyl-lignin. COMT downregulation also 
leads to the incorporation of 5'-hydroxy-guaiacyl units in 
lignin. However, syringyl-lignin does not exist in conifers, 
and we found no studies that show an effect of COMT 
downregulation on 5'-hydroxy-guaiacyl production in coni- 
fers. Because guaiacyl-lignin is the dominant lignin type in 



conifers, a decrease of COMT expression level could be as- 
sociated with a decrease of lignin synthesis. 

The upregulated putative laccase enzyme belongs to the 
wide super family of the multicopper oxidase (Nakamura 
and Go 2005). In plants, some laccase enzymes are involved 
in lignin biosynthesis, although they have a large spectrum 
of substrates and form a large family of genes. In loblolly 
pine, eight laccase genes have been described and two 
of them have been functionally characterized (Sato et al. 
2001; Sato and Whetten 2006). Both enzymes were able 
to oxidize coniferyl alcohol and produce dimers of coniferyl 
alcohol and as a consequence are involved in lignin 
biosynthesis. 

Two other upregulated genes in our constitutive samples 
are annotated as putative phenylcoumaran benzylic ether 
reductase. Phenylcoumaran benzylic ether reductases are 
involved in phenolic secondary metabolism and convert 
8'5'-linked lignin dihydroconiferyl alcohol into isodihydro- 
dehydrodiconiferyl alcohol by the reduction of benzylic 
ether functionality (Gang et al. 1999). A previous study 
showed that a phenylcoumaran benzylic ether reductase 
is involved in induced conifer defense following either me- 
chanical wounding or weevil attack (Lippert et al. 2007). 

The upregulated gene annotated as putative UDP-gluco- 
syltransferase plays an important role in lignin biosynthesis. 
After their biosynthesis, the monomers of lignin (i.e., p-cou- 
maryl, coniferyl, and sinapyl alcohols according to plant spe- 
cies) have to be translocated to the cell wall for the next 
oxidation step of lignin biosynthesis. The 4-0-(3-D-glucosides 
of cinnamyl alcohols have been considered as the transport 
forms of coniferyl and sinapyl alcohols. A UDPGxoniferyl al- 
cohol glucosyltransferase from Pinus strobus has been able 
to convert cinnamyl aldehydes as well as coniferyl and dihy- 
droconiferyl alcohols into their corresponding 0-(3-D-gluco- 
sides in vitro (Steeves et al. 2001). However, because 
coniferyl and sinapyl alcohols might be able to freely diffuse 
through the plasma membrane, it has been suggested that 
these glucosides play no role in monolignol export for de- 
velopmental lignin (Boija and Johansson 2006; Vanholme 
et al. 2008). Another noteworthy gene is annotated as pu- 
tative MYB16, a member of the family of transcription fac- 
tors. MYB16 belongs to the R2R3-MYB family and was 
shown to accumulate transiently in response to wounding 
in white spruce (Bedon et al. 2010) 

At least two genes are annotated within the flavonoid 
metabolism. First, an upregulated gene was annotated as 
a putative flavonoid 3'-monooxygenase that belongs to 
the cytochrome P450 superfamily. This gene is involved in 
central flavonoid metabolism, the leading precursors of fla- 
vones, anthocyanins, and proanthocyanidins pathways 
(Winkel-Shirley 2001). Anthocyanins can play various roles, 
including the resistance mechanisms toward insect pests 
(Steyn et al. 2002). The second gene within the flavonoid 
metabolism is downregulated and annotated as a putative 
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flavonol 4'-sulfotransferase. Ralph et al. (2006) found that 
several genes of flavonoid metabolism, including a Flavonoid 
3'-monooxygenase (= hydroxylase), are upregulated after 
white pine weevil herbivory, mechanical wounding, or 
western spruce budworm (Choristoneura occidentalis, 
Lepidoptera) feeding. 

Many Stress-Related Proteins Exist for Weevil 
Resistance 

Our study shows that 1 5 of the 26 putative sHSP and several 
other stress-related genes are downregulated in resistant 
trees. sHSP belong to a large family of proteins. They are 
highly variable, but they share a conserved oc-crystallin do- 
main of approximately 1 00 residues (Caspers et al. 1 995; de 
Jong et al. 1998; Fu et al. 2006). sHSP are classified into at 
least 1 1 subfamilies localized in different cell compartments: 
cytosol, mitochondria, chloroplasts, endoplasmic reticulum, 
and peroxisome (Vierling 1991b; Helm et al. 1993; Waters 
etal. 1996; Scharf etal. 2001; Siddique et al. 2003; Ma etal. 
2006; Waters et al. 2008). The 15 downregulated putative 
sHSP belong to class I, class II, class III, chloroplastic endo- 
plasmic reticulum, or cannot be assigned with confidence 
to a known class. The role of sHSP has been widely studied 
in plants. They are involved in plant response to various kinds 
of stress such as heat, cold, drought, heavy metals, salinity, 
oxidative, and osmotic stress (Vierling 1991a; Waters et al. 
1996; Wang et al. 2004; Haslbeck et al. 2005; Sun and 
MacRae 2005; Nakamoto and Vigh 2007). sHSP are also 
involved in normal development of plants, during embryo 
development, seed germination, somatic embryogenesis, 
pollen development, and fruit maturation (Sun et al. 
2002 and references therein). sHSP usually play a protection 
role (Haslbeck et al. 2005; Nakamoto and Vigh 2007). They 
can form stable complexes with denatured proteins to pre- 
vent its aggregation. sHSP also form soluble aggregates with 
substrate proteins, creating a transient reservoir of sub- 
strates. Release and refolding of both complexes and aggre- 
gates need the cooperation of ATP-dependent chaperone 
systems. sHSP also play a role in membrane quality control 
and are potential membrane stabilizing factors. 

Several sHSP were previously shown to be involved in co- 
nifer defense. Lippert et al. (2007) showed that weevil feed- 
ing induces the overexpression of seven sHSP at the protein 
level (up to 6-fold induction) in Sitka spruce. They also 
showed that transcript and protein expression levels are 
not correlated as six of the seven sHSP corresponding tran- 
scripts are not upregulated following weevil feeding. The 2- 
fold upregulation of the seventh sHSP transcript (class I) is 
comparable to the upregulation of the associated protein. 
Nevertheless, they observed that all the seven sHSP 
transcripts are constitutively expressed to high levels in bark 
tissue. Such constitutive expression of sHSP has also been 
observed in Ar. thaliana (Siddique et al. 2008), but the 



constitutive role of sHSP remains unknown. The results of 
Lippert et al. (2007) suggested that sHSP transcripts accumu- 
late in transient stocks and that sHSP expression is post- 
transcriptionally controlled. Recent studies have shown that 
RNA-binding proteins can regulate the stability, translation, 
or localization of mRNA (Glisovic et al. 2008; Hogan et al. 
2008; Babitzke et al. 2009). sHSP activity is also regulated at 
the protein level by phosphorylation or oligomer reorgani- 
zation. As a consequence, the expression levels of sHSP tran- 
scripts do not necessarily correlate with the sHSP expression 
at the protein level. sHSP may not play a role in constitutive 
defense and, in fact, may be involved in induced defense, 
among other biological processes. However, the test of this 
hypothesis needs a time-series comparison of both the 
transcriptome and the proteome after induction (e.g., 
weevil feeding), based on both susceptible and resistant 
strains of spruce. Together with 1 5 putative sHSP, 1 2 putative 
stress-related proteins are constitutively upregulated in 
susceptible trees. Their potential role is yet to be discovered. 

Phenotype Prediction and Efficiency of the 
Approach 

As in previous studies based on morphological features or 
terpene contents (Tomlin et al. 1997; Tomlin and Borden 
1997a; Alfaro et al. 2004), our goal was to determine if 
the transcriptome profiling is able to predict resistance levels 
in interior spruce. To determine whether the observed gene 
expression profiles corresponded to the observed pheno- 
type (i.e., resistant/susceptible), we performed a hierarchical 
clustering (fig. 4). Although the individuals clustered into 
two groupings, they did not match with the phenotype clas- 
sification. One cluster contained 11 susceptible trees and 
a second cluster contained the remaining trees, that is, both 
susceptible and resistant trees. The heat map clearly shows 
that 1 1 susceptible trees have a distinct profile of gene ex- 
pressions compared with the other 29 trees. Therefore, it 
might be possible to identify certain susceptible phenotypes 
by analyzing the transcriptome profiles, but it will not be 
possible to identify resistant trees with a high degree of cer- 
tainty using this approach. Four hypotheses could explain 
this pattern but at least three of them can be rejected. 

First, the resistance levels might be inaccurately assessed 
for some progenies. The family size of all the examined trees 
varied between 14 and 175 trees (see supplementary table 
S1, Supplementary Material online). Among the families 
used in the transcriptome comparison, six families (five sus- 
ceptible and on resistant) contained fewer than 80 individu- 
als: S-1 65-65, S-161-60, S-166-130, S-170-107, S-179-105, 
and R-1 1-19 (respectively, 42, 41, 30, 63, 14, and 42 trees). 
Four of them are considered susceptible and clustered with 
resistant trees in the cluster #2. Consequently, the assessment 
of the resistance levels of these progenies might be question- 
able. However, this does not explain why susceptible 
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progenies (with more than 80 individuals) S-1 76-1 33, S-1 56- 
103, S-1 58-1 31, S-1 67-95, and S-1 73-1 17 cluster with 
resistant trees. However, the original assessment of damage 
was based on natural levels of weevil attack. Attack patterns 
are rarely uniform in the wild, and all trees do not have 
the same probability of attack (He and Alfaro 1997). There- 
fore, some of the undamaged trees could have been 
"escapes" and never subject to attack, leading to some bias 
in the resistance levels assessment, particularly in the small 
progenies. 

Second, the differences in the observed damages caused 
by weevils can be explained by environmental factors such 
as growth conditions. This hypothesis seems improbable be- 
cause all the parent trees were collected within the same 
region (Prince Georges area) and the progenies were ran- 
domly mixed across several stands. All of them were grown 
in the same standard conditions. Moreover, as the trees used 
for gene expression profiling were grafted on the same root- 
stock, we do not expect high difference due to misadapta- 
tion to local soil conditions. 

Third, as the collected seeds were open pollinated in the 
wild, we know only the mother and have no information 
about the fathers of the progenies used for resistance scor- 
ing. This may induce a bias if parents have very different lev- 
els of resistance. However, a previous study has shown 
a high family heritability {h 2 = 0.70) in a similar experiment 
design (King et al. 1997), and crosses between susceptible 
and resistant trees would lead to intermediate levels of re- 
sistance (Alfaro et al. 2004). As a consequence, a bias in- 
duced by the uncertainty of fatherhood of the progenies 
seems improbable. 

Finally, the resistance or susceptibility may be based on 
several different strategies, involving different sets of genes. 
In this case, our experimental design does not allow us to 
identify genes involved only in rare strategies. If resistance 
can be associated with, for example, ten different profiles of 
gene expression, we can expect only a few trees for each 
strategy to be present in our sampling. In such a case, 
the differences in gene expression profiles will be confused 
with individual variations because we did not classify 
the trees according to their strategy but according to their 
phenotype. 

Supplementary Material 

Supplementary table S1 and S2 are available at Genome Biology 
and Evolution online (http://www.gbe.oxfordjournals.org/). 
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